arXiv: 1505.02225vl [physics.geo-ph] 9 May 2015 


Aftershocks and Omori’s law in a modified Carlson-Langer model with nonlinear 

visco-elasticity 

Hidetsugu Sakaguchi and Kazuki Okamura 
Department of Applied Science for Electronics and Materials, 

Interdisciplinary Graduate School of Engineering Sciences, 

Kyushu University, Kasuga, Fukuoka 816-8580, Japan 

A modified Carlson-Langer model for earthquakes is proposed, which includes nonlinear visco¬ 
elasticity. Several aftershocks are generated after the main shock owing to the damping of the 
additional visco-elastic force. Both the Gutenberg-Richter law and Omori’s law are reproduced in 
a numerical simulation of the modified Carlson-Langer model on a critical percolation cluster of a 
square lattice. 
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I. INTRODUCTION 


Earthquakes follow various types of power laws. Gutenberg and Richter found that the slip-size of faults or the 
seismic moment S obeys a power law: [l| 

P{M) - 10“'’^, (1) 

where the magnitude M is defined as M = (2/3) logj^Q S. The exponent b takes a value ranging from 0.8 to 1.2, and 
the most typical value is 5 = 1. After a large main shock, many aftershocks occur. Omori found that the frequency 
of aftershocks decays in a power law: [3, 
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where r is the elapsed time from the main shock, and the exponent p is expected to take a value between 0.9 and 1.4. 
In the original Omori’s law, p is set to be 1. Omori’s law implies that the frequency of aftershocks decays slowly or 
aftershocks can occur even after a rather long time following the main shock. The mechanisms of earthquakes and 
aftershocks have been studied by many authors, but are not completely understood. There are several deterministic 
models for earthquakes. The power laws of earthquakes might be understood from the chaotic dynamics of the 
nonlinear equations. Burridge and Knopoff proposed a block-spring model for earthquakes Later, Carlson and 
Langer used a type of velocity-weakening friction and found a power law in the small magnitude range Many 
authors studied intensively the Carlson-Langer model and its modification. Vieira, Vasconcelos, Nagel pointed out a 
transition in the magnitude distribution by changing a parameter of the velocity-weakening friction [6| . On the other 
hand, there are various models in which the power laws originate from the fractal geometry of faults or asperities 
of plate boundaries. The size distribution of earthquakes has been explained by several authors using such fractal 
geometry Q. We studied the block-spring system on critical percolation clusters with fractal geometry, and showed 
that the exponent of the power law is about 1, which is close to the exponent in the original Gutenberg-Richter law Q. 
There are also several stochastic cellular automaton models for earthquakes. For example, Bak and Tang proposed a 
cellular automaton model for earthquakes which exhibits self-organized criticality Q. Olami, Feder, and Christensen 
proposed a modified cellular automaton model with some dissipation 0- 

The difference between main shocks and aftershocks is not clear in the Carlson-Langer model or the Bak-Tang 
model. Ito and Matsuzaki proposed a modified model of Bak-Tang model and rmroduced Omori’s law for aftershocks 
by randomly perturbing certain regions that slipped during the main shock [ll|. Dieterich studied time-dependent 
friction in rocks and suggested a relationship between the visco-elasticity and Omori’s law [T^ . Hainzl et al. included 
some visco-elasticity in the Olami-Feder-Christensen model and reproduced a version of Omori’s la w 11^ . Jagla et 
al. studied a depinning model of viscoelastic interfaces and reproduced a type of Omori’s law [l3 . Il5j| . However, 
the physical meaning of models constructed by Ito-Matsuzaki and Hainzl et al. remains unclear. Jagla’s depinning 
model includes randomness and the relationship between the depinning model and earthquakes is not obvious. In this 
paper, we study a modified version of the Carlson-Langer model with visco-elastic elements, and discuss Omori’s law 
of aftershocks. Our model is based on mechanics or Newton’s equation of motion, and does not include randomness or 
external noise in contrast to the models by Hainzl et al. and Jagla et al. The time evolution is continuous during slip 
events, but stick intervals between successive slips are theoretically evaluated and numerical simulation for the stick 
intervals is skipped. A rather long-time numerical simulation becomes possible owing to the reduction in computation 
time. 



FIG. 1: Schematic figure of a modified Carlson-Langer model with visco-elasticity. 


II. THE CARLSON-LANGER MODEL WITH VISCOELASTICITY 


The Carlson-Langer model is a block-spring model. In a one-dimensional model, N blocks of mass to = 1 are 
coupled with the nearest neighbor blocks using a spring with a spring constant k. In the original model, each block is 
further coupled with a spring to a rigid plate moving with a constant velocity vq. The model equation is expressed as 

(Pxi 

= k{xi+i - 2xi + Xi-i) -I- {vot - Xi) - (3) 

where Vi denotes the velocity of the block dxi/dt^ k is the spring constant between neighboring blocks, and 4‘i'^i) 
represents the kinetic friction. The maximum static friction is set to be 1. In the Carlson-Langer model, the kinetic 
friction is assumed to be 
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for u > 0. Two parameters a and a characterize the velocity-weakening friction. The parameter a represents the 
difference between the maximum static friction 1 and the kinetic friction at the limit of n = 0. The parameters k 
and (T are fixed to be fc = 16 and a = 0.01 in our numerical simulations for the sake of simplicity. The parameter 
a expresses the velocity-weakening rate of the kinetic friction, i.e., d(j)(y)/dv = —2a at the limit of u = 0. When Vi 
decreases to 0, the slip state changes to the stick state. The stick state is assumed to be maintained until the applied 
force reaches the maximum static friction 1. It is assumed that the reverse motion does not occur Q. 

In the modified model, a dashpot and spring are added in parallel to the spring coupled with the rigid plate pulled 
by velocity vq as shown in Fig. 1. An additional force /i < 0 is applied to the ith block by the additional dashpot 
and spring after the block starts to slip. The additional force works as a braking force to the slipping block. Various 
models including viscoelasticity are possible. For example, the dashpot and the spring can be added in parallel to 
the spring connecting neighboring blocks. These models might be interesting but we investigate only the simplest 
model shown in Fig. I in this paper. One reason is that the role of the viscoelasticity as a brake force is rather clear 
and anther reason is that the dynamics in the stick phase is simple and analytically calculated as explained in later 
paragraphs. 

The force satisfies fi = Ky^ where —yi represents the contraction of the additional spring. For a linear dashpot, 
the force fi satisfies fi = —f]{—Ui) where Ui < 0 denotes the slip velocity at the dashpot. Because Omori’s law is 
not reproduced in a model using the linear dashpot as shown later, we use a more general model where fi satisfies 
fi = —r]{—Ui)'^. A power-law relation is assumed between the slip velocity and resistance force, and 7 = 1 corresponds 
to the linear dashpot. A power-law model for the relationship between shear velocity and shear stress is one of the 
typical models for non-Newton fluids. The power-law model of 0 < 7 < I is applied to non-Newton fluids such as 
polymeric fluids, tooth paste and chocolate. The slip velocity ut is expressed as Ui = —{—fi/yY, where /3 = I/ 7 . The 
velocity difference between the upper plate and the Ah block is written as vg — dxi/dt which is equal to dyi/dt + Ui. 
Then, Xi and fi satisfy 
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We further assume that Vq is infinitesimally small, and Vg is set to be zero during the slip phase. In the stationary 
stick phase, fi = 0, and Vi = 0. First, we calculate the maximum value F oi Fi = k{xi+i — 2xi -\- Xi-i) + (vgt — Xi) 













3 


among the N blocks. If the ith block takes the maximum value, the time is advanced from t to t + (1 — F)/vo. The 
pulling force at the ith block exceeds the maximum static friction 1, and the ith block starts to slip. The neighboring 
blocks can slip together owing to the elastic coupling of spring constant k. The time evolution of the slip process 
obeys Eq. (5). Because we assume vq is infinitesimally small, two slip events do not occur simultaneously in distant 
places. When the first slip event (main shock) starts, dxi/dt > 0 for the block i, and fi becomes negative. That 
is, the braking force fi works for the slipping block. After the main shock ends and the stick phase starts, the slip 
velocity dxi/dt becomes zero, then /i < 0 increases toward zero by the last term K{—fi/ri)^ > 0 owing to the effect 
of the dashpot. That is, the braking force dampens and the total pulling force to the block tends to increase over 
time. It is possible that some blocks start to slip again after a certain time following the main shock. This is an 
aftershock. Aftershocks can occur repeatedly. After the sequence of aftershocks, all ffs decay to zero. The period 
from the start of a main shock and the end of the sequence of aftershocks is one shock event. At the end of the shock 
event, Fi = k{xi+i — 2xi + Xi-i) + (vot — Xi) is smaller than 1 and fi — 0 for all i’s. Then, we again calculate the 
maximum value F oi Fi = k{xi+i — 2xi + Xi-i) + {v^t — Xi) among the N blocks, and the upper plate is slowly moved 
until another block begins to slip. The time interval At until the starting time of the next main shock is calculated 
by foAt = \ — F. Similarly, we can calculate the start time of aftershocks in the sequence of aftershocks. In the stick 
phase after an aftershock, fi{t) can be explicitly expressed as 


+ (1 - (6) 

■qp 

for /3 ^ 1 because uq = 0 and dxi/dt = 0 after the final time to of the aftershock. From Eq. ( 6 ), a possible start time 
ti of the next aftershock for the ith block is evaluated as 


ti — to + 
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because Fi + fi (t) = 1 at the start time ti of the next aftershock. The minimum value tg of ti is the actual start time 
of the next aftershock, and fi{ts) at t = tg is evaluated using Eq. ( 6 ). If all Ffs are smaller than 1, the sequence of 
aftershocks ends. In the case of the linear visco-elasticity: 7 = /3 = 1, Eqs.( 6 ) and (7) are replaced by 
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The equation of motion Eq. (5) with Eq. (4) is used only in the slip phases, and the time intervals in the stick phases 
between main shocks or between two aftershocks are theoretically calculated using Eq. (7) or Eq. (9). We performed a 
numerical simulation repeating these processes. A characteristic of our model is that a sequence of aftershocks occurs 
after each main shock and the difference between main shocks and aftershocks is clear. Because vq is assumed to 
be infinitesimally small, the time of a sequence of aftershocks is infinitesimally small compared with the interval of 
successive main shocks, which is an order of I/vq. 


III. NUMERICAL RESULTS IN ONE DIMENSION 

First, we show the numerical results of the modified Carlson-Langer model of ten blocks to understand the role 
of aftershocks. Figure 2(a) shows the time evolution of the sum S' = ^ of slip distances at 7 = 1/2 (/3 = 2), 
a = 1.5, AT = 4, and 6 = Kjrf = 1. In the time scale of O(l/uo), aftershocks occur at the same time as the main 
shock. The vertical lines show aftershocks. The largest value of S at each slip time corresponds to the main shock. 
A slip event without aftershocks is shown as a rhombus without vertical lines. The number of rhombi at each slip 
time is the sum of one (main shock) plus the number of aftershocks. The number of aftershocks takes a various value 
ranging from 0 to 6 . Figure 2(b) shows time evolution of the profile of Xi for the ten-block system using the same 
parameters. The light solid lines denote the profiles of Xi after the main shocks and the dark dashed lines denote 
the profiles of Xi after the aftershocks. Figure 2(c) is the enlargement of Fig. 2(b) for 424.26 < x < 424.272. After 
a large main shock, a dent structure is created, which is composed of delayed blocks. Similar dent structures were 
also found in the original Carlson-Langer model Q. Pulling forces are relatively strong in these dent regions, because 
Xi+i — 2xi + Xi-i is positive and large. As the braking forces fi weaken, aftershocks tend to occur near these dent 
structures. Aftershocks make the profile of Xi flat. When the profile of Xi becomes sufficiently flat, a simultaneous 





4 


(a) 



424.45 

424.4 

424.35 

424.3 

424.25 

0123456789 

i 


(0 



i 


(b) 


“I—I—I I I I I I 


FIG. 2: (color online) (a) Time evolution of the sum S of slips for ten-block systems at 7 = 1/2 (/3 = 2), a = 1.5, K = 4, and 
6 = K/r}^ = 1. (b) Time evolution of the profile of Xi for ten-block systems using the same parameters. The light solid lines 
denote the profiles of Xi after main shocks and the dark dashed lines denote the profiles of Xi after aftershocks, (c) Enlargement 
of (b) for 424.26 <x< 424.272. 
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FIG. 3: (a) Semi-logarithmic plot of the probability distribution of the elapsed time r of aftershocks at 7 = l,Q = l,E' = 0.5, and 
S — K/rj = 0.05. (b) Double-logarithmic plot of the probability distribution of r at 7 = 1/2, a = 1, K = 4, and <5 = K/tf — 1. 
(c) Double-logarithmic plot of the probability distribution of r at 7 = 1/3, q = 0.5, if = 0.4, and 5 = K/rj^ = 10000. (d) 
Double-logarithmic plot of the probability distribution of r at 7 = 1/4, a = 0.5, K = 0.4, and 5 = Kjrj'^ = 100000. 


slip or a large shock occurs next. Through chaotic dynamics in the simultaneous slip event, the spatial fluctuation 
increases over time and another dent structure appears after the large-scale shock. Similar intermittent dynamics 
were also observed in the original Carlson-Langer model In the original Carlson-Langer model, small slip evens 
(small-scale main shocks) frequently occurs near the dent structures. In our modified model, aftershocks play a role 
in these small slip events in the original Carlson Danger model. 

We investigated the statistical properties of aftershocks in a one-dimensional system of = 1000. We performed 
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FIG. 4: Probability distribution of the magnitude at (a) a = 0.3, (b) 0.4, and (c) 0.5 for 7 = 1/4, K — 0.4, and 5 = K/tf = 
100000. The dashed line in Fig.4(b) denotes a curve of P{M) oc 10“^’^“^. 


numerical simulations by changing parameters 7, a,i4r and 5 = . We present typical examples in this paper. 

Figure 3(a) shows a semi-logarithmic plot of the probability distribution of the time t when aftershocks occur after the 
main shocks at 7 = 1, a = 1, iF = 0.5, and 5 = iF/ry = 0.05. The dashed line is n(r) = 0.048 exp(—0.048T), although 
the dashed line overlaps with the solid curve and is hardly visible. In this case of linear visco-elasticity 7 = 1 , 
the probability distribution obeys an exponential distribution. This implies that the simple linear model does not 
reproduce Omori’s law of aftershocks. Figure 3(b) shows a double-logarithmic plot of the probability distribution of r 
at 7 = 1/2, a = 1, iF = 4, and 8 = K/rf = 1. The dashed line is n(T) = 356/(373 -I- r)^. Figure 3(c) shows a double- 
logarithmic plot of the probability distribution of r at 7 = 1/3, a = 0.5, K = 0.4 and <5 = K/rf = 10000. The dashed 
line is n(T) = 0.98/(31 -I- r)^'^. Figure 3(d) shows a double-logarithmic plot of n(r) at 7 = 1/4, a = 0.5, iF = 0.4, 
and S = K/rj* = 100000. The dashed line is n(T) = 0.18/(310 -I- t)^'°^. Similarly, the probability distribution at 
7 = 3 / 4,0 = 1, K = 1,S = Kfr]^/'^ = 0.1 is also calculated and the approximate function n{T) = 60000/(350-I-r)^'® 
is obtained. The probability distributions approximately obey Omori’s law in the nonlinear visco-elastic models. The 
exponents p are evaluated at p = 2.8, 2,1.3 and 1.02 for 7 = 3/4,1/2,1/3 and 1/4, respectively. The exponent p is 
interpreted as p = 00 for the linear case of 7 = 1. These results suggest that the exponent p decreases with 7. For 
7 = 1/4, the exponent p is close to the exponent p = 1 of original Omori’s law. 

Figures 4 shows the probability distribution P{M) of the magnitude at (a) a = 0.3, (b) 0.4, and (c) 0.5 for 
7 = 1/4, iF = 0.4, and S = K/rf = 100000. Figure 4(c) is a magnitude distribution for the numerical simulation 
shown in Fig. 3(d). In the original Carlson-Langer model, the magnitude distribution has a hump structure around 
a rather large value of M at a > 1 and decays rapidly for large M at a < 1. The magnitude distibution exhibits a 
power law behavior near a = 1. Even in our model with nonlinear viscoelastisity, the magnitude distribution decays 
rapidly for large M at small a as shown in Fig. 4(a) and has a long tail at large a as shown in Fig.4(c), although 
there is a shoulderlike structure near M = —0.7 in the magnitude distributions. The magnitude distribution exhibits 
a power-law-like behavior at a = 0.4 for M > —0.5, although the linearity is not so good as compared to the models 
on a percolation cluster as shown in Fig. 6(a) and Fig. 7(a). The exponent b is evaluated as 1.2, which is close to the 
b value of the Gutenberg-Richter law for earthquakes. 


IV. NUMERICAL RESULTS ON A PERCOLATION CLUSTER 

Next, we investigated the statistical properties of aftershocks on a critical percolation cluster on a square lattice 
of 150 X 150. The critical percolation cluster was constructed by the site percolation on a square lattice. Only the 
neighboring occupied sites are connected by springs of a spring constant k = 16. We assumed that blocks are located 
only in the largest mutually-connected percolation cluster. In the previous papen we performed a numerical simulation 
of the original Carlson-Langer model on the same critical percolation cluster Q. We show numerical results for the 
modified Carlson-Langer model on the same percolation cluster at 7 = 1/2, a = 1.2, K = 1, and S = K/rf = 0.05. 
Figure 5(a) shows the time evolution of the total slip size S of the main shocks. Figure 5(b) shows the time evolution 
of the aftershocks after a main shock at v^t = 3.0812. The frequency of aftershocks decreases with time. Figure 5(c) 
shows epicenters of the main shock (rhombus) and aftershocks (pluses). The epicenters of aftershocks locate around 
the epicenter of the main shock. In real earthquakes, it is observed that epicenters of aftershocks is locally distributed 
around the epicenter of the main shock. Here, the epicenters are defined as the position of the block which slips first 
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FIG. 5: (color online) (a) Time evolution of main shocks at 7 = 1/2,a = 1.2, K — 1, and 5 — K/rf = 0.05. (b) Time evolution 
of aftershocks after the main shock at vqI = 3.0812. (c) Distribution of the epicenters of a main shock (rhombus) and the 
aftershocks (pluses). The critical percolation cluster is expressed using dots, (d) Relationship between the slip size S of the 
main shock and the number Na of aftershocks after each main shock. The dashed line is Na = 


(a) (b) (c) 



FIG. 6 : (a) Probability distribution of the magnitude M = (2/3) logos' at 7 = 1/2, a = 1.2, K — 1, and 5 — K/rf = 0.05. 
The dashed line denotes an exponential distribution P(M) oc 10”^’^ . (b) Probability distribution of r. The dashed line is 

n(r) = 13.3/(680-|-r)^ ®. (c) Probability distribution of r for M > 0 and M < 0. The dashed lines are n(r) = 7.4/(2560 + r)^'®® 
and n(T) = 598/(916 -I- 


for each slip event. Figure 5(d) shows a relationship between the total slip size S of the main shock and the number 
Na of aftershocks after each main shock. The number of aftershocks increases with the magnitude of the main shock. 
The dashed line is Na = 95'° ®^. 

Figure 6 (a) is the probability distribution of the magnitude M = (2/3)logjQS' at 7 = 1/2, a = 1.2, K = 1, and 
S = K/rj^ = 0.05. The dashed line denotes a distribution P{M) cx 10“^'^^. A good agreement is seen for M < 1.5. A 
type of Gutenberg-Richter law is satisfied. Figure 6 (b) shows the probability distribution of the elapsed time r of the 
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FIG. 7: (a) Probability distribution of the magnitude M — (2/3) logjn S' at 7 = 1/4, a = 0.95, K — 0.1, and 5 = Kjvi^ — 100000. 
The dashed line denotes an exponential distribution P{M) oc 10”^' (b) Probability distribution of r. The dashed line is 

n(r) = 0.09/(48 +t)^ °^. (c) Probability distribution of r for M > 0 and M < 0. The dashed lines are n{T) = 0.0068/(112+r)^'^ 
and n(r) = 0.0068/(215 + . 


aftershocks after the main shocks. The dashed line shows n(r) = 13.3/(680 + r)^-®. Omori’s law of p = 1.5 is obtained. 
We studied the difference of the the probability distributions of r after the main shocks of larger magnitude M > 0 
and smaller magnitude M < 0. Figure 6 (c) shows the two probability distributions. Both distributions obey Omori’s 
law but the exponent p is 1.35 for M > 0 and 1.95 for M < 0. We also performed a numerical simulation on the same 
percolation cluster on a square lattice of 150 x 150 at 7 = 1/4, a = 0.95, iF = 0.1, and 6 = = 100000. Figure 

7(a) is the probability distribution of the magnitude M = (2/3) log^Q 5”. The dashed line denotes an exponential 
distribution P{M) oc 10“^ °®^. A good agreement is observed for M < 1.5. The distribution P{M) oc lO”**^ with 
6 = 1 corresponds to the original Gutenberg-Richter law. Figure 7(b) shows the probability distribution of r. The 
dashed line shows n(r) = 0.09/(48 -I- r)^ °^. Omori’s law of p = 1.02 is obtained. The Gutenberg-Richter law and 
Omori’s law close to the original versions are obtained at this parameter set. Figure 7(c) shows the two probability 
distributions of r for M >0 and M < 0. Both distributions obey Omori’s law but the exponent p is 1.1 for M > 0 and 
0.75 for M < 0. We do not understand the reason but the exponent seems to depend on the magnitude of the main 
shock and the exponent p is larger for larger M in our model. We also performed a similar numerical simulation for a 
linear visco-elastic model of 7 =lata = l,Ar = 0.5, and S = K/rj = 1. We found that the probability distribution of 
r obeys an exponential function in the linear model, which is the same as the behavior of the one-dimensional model 
shown in Fig. 3(a). These numerical results suggest that the exponent p of Omori’s law does not strongly depend on 
the dimension. 


V. SUMMARY 

We have proposed a modified Garlson-Langer model which can generate aftershocks, and performed numerical 
simulations. The visco-elastic element works as a brake in the slip process of blocks. Owing to the damping of 
the braking force, the forward tensile force for blocks increases, and aftershocks are induced. By using nonlinear 
visco-elastic elements, Omori’s law has been reproduced. On the other hand, Omori’s law could not be reproduced 
in simple linear visco-elastic models. Both the Gutenberg-Richter law with exponent close to 1 and Omori’s law with 
exponent close to 1 are observed in the numerical simulation of the modified Garlson-Langer model with 7 = 1/4 on 
a critical percolation cluster on a square lattice. 
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